{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import matplotlib.pyplot as plt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!ls -sht N_40*.dat"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "134M N_40__T1_4__T2_4__dm_0.44.dat\r\n",
        "133M N_40__T1_4__T2_12__dm_0.44.dat\r\n",
        "133M N_40__T1_4__T2_8__dm_0.44.dat\r\n",
        "134M N_40__T1_4__T2_4__dm_0.22.dat\r\n",
        "133M N_40__T1_4__T2_8__dm_0.22.dat\r\n",
        "133M N_40__T1_4__T2_12__dm_0.22.dat\r\n",
        "133M N_40__T1_4__T2_8__dm_0.11.dat\r\n",
        "133M N_40__T1_4__T2_12__dm_0.11.dat\r\n",
        "134M N_40__T1_4__T2_4__dm_0.11.dat\r\n",
        "133M N_40__T1_4__T2_8__dm_0.0078.dat\r\n",
        "132M N_40__T1_4__T2_12__dm_0.0078.dat\r\n",
        "134M N_40__T1_4__T2_4__dm_0.0078.dat\r\n",
        "134M N_40__T1_4__T2_4__dm_0.0.dat\r\n",
        "364K N_40__T1_4__T2_4__dm_0.2.dat\r\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "arquivos=!ls -t N_40__T1_4__T2_8__dm_*.dat"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "arquivos"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 3,
       "text": [
        "['N_40__T1_4__T2_8__dm_0.0.dat',\n",
        " 'N_40__T1_4__T2_8__dm_0.0078.dat',\n",
        " 'N_40__T1_4__T2_8__dm_0.11.dat',\n",
        " 'N_40__T1_4__T2_8__dm_0.22.dat',\n",
        " 'N_40__T1_4__T2_8__dm_0.44.dat']"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dm=[0.0, 0.0078, 0.11, 0.22, 0.44]\n",
      "#dm=[0.11, 0.22, 0.44]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n",
      "font = {'family' : 'normal',\n",
      "         'weight' : 'bold',\n",
      "         'size'   : 50}\n",
      "plt.rc('font', **font)\n",
      "plt.rc('text', usetex=True)\n",
      "\n",
      "transient= 100000\n",
      "m_every = 5000\n",
      "dm=[0.0, 0.0078, 0.11, 0.22, 0.44] #range(4,48,4)\n",
      "#dm=[0.11, 0.22, 0.44]\n",
      "Npart=41\n",
      "arquivos = [[\"N_\"+str(Npart)+\"__T1_4__T2_8__dm_\"+str(i)+\".dat\", i] for i in dm]\n",
      "Temp=np.zeros(len(arquivos)*Npart).reshape(len(arquivos),Npart)\n",
      "for w in range(len(arquivos)):\n",
      "  data = np.loadtxt(arquivos[w][0])\n",
      "  ttime=np.array(data[:,0])\n",
      "  vel=np.array([data[:,2*i+1] for i in range(Npart)])\n",
      "  #pos=np.array([data[:,2*i+2] for i in range(40)])\n",
      "  del data\n",
      "\n",
      "  Temp[w]=vel[:,transient:].var(axis=1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Correct temperature\n",
      "TempCorrct=np.zeros(len(Temp)*len(Temp[0])).reshape(len(Temp),len(Temp[0]))\n",
      "#Temp[0]*(dm[0])\n",
      "for w in range(len(Temp)):\n",
      "    for i in range(len(Temp[w])):\n",
      "        TempCorrct[w][i]=Temp[w][i]*(1+(i%2)*dm[w])\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for w in range(len(arquivos)):\n",
      "    plt.plot(TempCorrct[w],'o-', label=\"$\\delta m=\"+str(arquivos[w][1])+\"$\")\n",
      "#plt.plot(vel[:,transient:].mean(axis=1),'x-')\n",
      "plt.plot(XxX,YyY)\n",
      "plt.ylim([0,9])\n",
      "plt.xlabel(\"$n$\")\n",
      "plt.ylabel(\"$T$\")\n",
      "plt.grid()\n",
      "plt.legend(loc=8, prop={'size':20}, mode=\"expand\", ncol=6)\n",
      "plt.subplots_adjust(bottom=0.15)\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def Tk(x,T1,T2,Np):\n",
      "    return (T1**(3.0/2.0)-(T1**(3.0/2.0)-T2**(3.0/2.0))*x/Np)**(2.0/3.0)\n",
      "XxX=range(41)\n",
      "YyY=range(41)\n",
      "for i in range(41):\n",
      "    YyY[i]=Tk(XxX[i],2,8,41)\n",
      "plt.plot(XxX,YyY)\n",
      "#plt.plot(TempCorrct[4])\n",
      "plt.grid()\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.plot(TempCorrct[0])\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "heading",
     "level": 6,
     "metadata": {},
     "source": [
      "Analisis do transiente para $\\delta m$ grande"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "arquivos[2][0]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 12,
       "text": [
        "'N_41__T1_4__T2_4__dm_0.44.dat'"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "data = np.loadtxt(arquivos[2][0])\n",
      "ttime=np.array(data[:,0])\n",
      "vel=np.array([data[:,2*i+1] for i in range(40)])\n",
      "pos=np.array([data[:,2*i+2] for i in range(40)])\n",
      "del data"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 6,
     "metadata": {},
     "source": [
      "Temporal Evolution"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig, axes = plt.subplots(2, sharex=True)\n",
      "transient= 50000\n",
      "m_every = 1000\n",
      "for i in range(40):\n",
      "    axes[0].plot(ttime[:transient], vel[i][:transient], markevery=m_every)\n",
      "    axes[1].plot(ttime[:transient], pos[i][:transient], markevery=m_every)\n",
      "# Make it beuty\n",
      "axes[0].set_ylabel(\"$v_x$\")\n",
      "axes[1].set_ylabel(\"$x$\")\n",
      "axes[1].set_xlabel(\"$t$\")\n",
      "axes[0].locator_params(axis='y',nbins=4)\n",
      "axes[1].locator_params(axis='y',nbins=4)\n",
      "axes[0].grid()\n",
      "axes[1].grid()\n",
      "plt.subplots_adjust(bottom=0.15)\n",
      "#plt.tight_layout()\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "AttributeError",
       "evalue": "'NoneType' object has no attribute 'tk'",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
        "\u001b[0;32m<ipython-input-14-a9bac6ba9bdd>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     15\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msubplots_adjust\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbottom\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.15\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     16\u001b[0m \u001b[0;31m#plt.tight_layout()\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 17\u001b[0;31m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshow\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
        "\u001b[0;32m/usr/lib/pymodules/python2.7/matplotlib/pyplot.pyc\u001b[0m in \u001b[0;36mshow\u001b[0;34m(*args, **kw)\u001b[0m\n\u001b[1;32m    137\u001b[0m     \"\"\"\n\u001b[1;32m    138\u001b[0m     \u001b[0;32mglobal\u001b[0m \u001b[0m_show\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 139\u001b[0;31m     \u001b[0m_show\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkw\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    140\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    141\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;32m/usr/lib/pymodules/python2.7/matplotlib/backend_bases.pyc\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, block)\u001b[0m\n\u001b[1;32m    107\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    108\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mis_interactive\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 109\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmainloop\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    110\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    111\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mmainloop\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;32m/usr/lib/pymodules/python2.7/matplotlib/backends/backend_tkagg.pyc\u001b[0m in \u001b[0;36mmainloop\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m     67\u001b[0m \u001b[0;32mclass\u001b[0m \u001b[0mShow\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mShowBase\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     68\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mmainloop\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 69\u001b[0;31m         \u001b[0mTk\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmainloop\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     70\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     71\u001b[0m \u001b[0mshow\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mShow\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;32m/usr/lib/python2.7/lib-tk/Tkinter.py\u001b[0m in \u001b[0;36mmainloop\u001b[0;34m(n)\u001b[0m\n\u001b[1;32m    326\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mmainloop\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    327\u001b[0m     \u001b[0;34m\"\"\"Run the main loop of Tcl.\"\"\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 328\u001b[0;31m     \u001b[0m_default_root\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtk\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmainloop\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    329\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    330\u001b[0m \u001b[0mgetint\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mint\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;31mAttributeError\u001b[0m: 'NoneType' object has no attribute 'tk'"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "heading",
     "level": 6,
     "metadata": {},
     "source": [
      "Impar number of particles"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "data = np.loadtxt(\"N_41__T1_4__T2_4__dm_0.44.dat\")\n",
      "ttime=np.array(data[:,0])\n",
      "vel=np.array([data[:,2*i+1] for i in range(40)])\n",
      "pos=np.array([data[:,2*i+2] for i in range(40)])\n",
      "del data"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "transient= 100000\n",
      "m_every = 5000\n",
      "Temp=vel[:,transient:].var(axis=1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "font = {'family' : 'normal',\n",
      "         'weight' : 'bold',\n",
      "         'size'   : 50}\n",
      "plt.rc('font', **font)\n",
      "plt.rc('text', usetex=True)\n",
      "\n",
      "transient= 100000\n",
      "m_every = 5000\n",
      "\n",
      "plt.plot(Temp,'o-', label=\"$\\delta m=\"+\"0.0\"+\"$\")\n",
      "#plt.plot(vel[:,transient:].mean(axis=1),'x-')\n",
      "#plt.ylim([0,17])\n",
      "plt.xlabel(\"$n$\")\n",
      "plt.ylabel(\"$T$\")\n",
      "plt.grid()\n",
      "plt.legend(loc=8, prop={'size':20}, mode=\"expand\", ncol=6)\n",
      "plt.subplots_adjust(bottom=0.15)\n",
      "plt.show()\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "/usr/lib/pymodules/python2.7/matplotlib/offsetbox.py:75: RuntimeWarning: divide by zero encountered in double_scalars\n",
        "  sep = (total - sum(w_list))/(len(w_list)-1.)\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "vel"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "array([[  2.00006000e+00,   2.00907000e+00,   2.02766000e+00, ...,\n",
        "          1.60043000e+00,   1.60043000e+00,   1.60043000e+00],\n",
        "       [ -1.63740000e-07,  -2.20844000e-03,  -1.10712000e-02, ...,\n",
        "         -1.17288000e+00,  -1.16568000e+00,  -1.15766000e+00],\n",
        "       [  0.00000000e+00,   3.93650000e-07,   7.31557000e-06, ...,\n",
        "         -1.43893000e+00,  -1.43818000e+00,  -1.44547000e+00],\n",
        "       ..., \n",
        "       [  0.00000000e+00,   2.07714000e-11,   1.46203000e-09, ...,\n",
        "         -6.39756000e-01,  -6.30591000e-01,  -6.22578000e-01],\n",
        "       [  0.00000000e+00,  -3.93650000e-07,  -7.31557000e-06, ...,\n",
        "         -2.54965000e-01,  -2.72087000e-01,  -2.87300000e-01],\n",
        "       [  1.63740000e-07,   2.20844000e-03,   1.10712000e-02, ...,\n",
        "          2.37272000e-01,   2.39041000e-01,   2.42366000e-01]])"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "heading",
     "level": 6,
     "metadata": {},
     "source": [
      "Compara\u00e7\u00e3o teorica"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ddd=\"N_40__T1_4__T2_8__dm_0.22.dat\"\n",
      "data = np.loadtxt(\"N_40__T1_4__T2_8__dm_0.22.dat\")\n",
      "ttime=np.array(data[:,0])\n",
      "vel=np.array([data[:,2*i+1] for i in range(40)])\n",
      "pos=np.array([data[:,2*i+2] for i in range(40)])\n",
      "del data"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "transient=100000\n",
      "Temp=np.zeros(40)\n",
      "Temp=vel[:,transient:].var(axis=1)\n",
      "TempCorrct=np.zeros(40)\n",
      "for i in range(40):\n",
      "    TempCorrct[i]=Temp[i]*(1+(i%2)*0.22)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def Tk(x,T1,T2,Np):\n",
      "    return (T1**(3.0/2.0)-(T1**(3.0/2.0)-T2**(3.0/2.0))*x/Np)**(2.0/3.0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def Tk(x,T1,T2,Np):\n",
      "    return (T1**(3.0/2.0)-(T1**(3.0/2.0)-T2**(3.0/2.0))*x/Np)**(2.0/3.0)\n",
      "\n",
      "XxX=range(40)\n",
      "YyY=range(40)\n",
      "for i in range(40):\n",
      "    YyY[i]=Tk(XxX[i],4,8,40)\n",
      "plt.plot(XxX,YyY)\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "font = {'family' : 'normal',\n",
      "         'weight' : 'bold',\n",
      "         'size'   : 50}\n",
      "plt.rc('font', **font)\n",
      "plt.rc('text', usetex=True)\n",
      "\n",
      "transient= 100000\n",
      "m_every = 5000\n",
      "\n",
      "plt.plot(TempCorrct[0],'o-', label=\"$\\delta m=\"+\"0.0\"+\"$\")\n",
      "plt.plot(XxX,YyY)\n",
      "#plt.plot(vel[:,transient:].mean(axis=1),'x-')\n",
      "#plt.ylim([0,17])\n",
      "plt.xlabel(\"$n$\")\n",
      "plt.ylabel(\"$T$\")\n",
      "plt.grid()\n",
      "plt.legend(loc=8, prop={'size':20}, mode=\"expand\", ncol=6)\n",
      "plt.subplots_adjust(bottom=0.15)\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}